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’  ABSTRACT 

'\ 

t.'' 

A  method  of  phase  comparison  time  delay  estimation  using  large 
time-bandwidth  product  signals  is  presented.  This  method  compares  the 
phases  of  the  matched  filters  for  each  channel,  and  it  is  shown  that  for 
signals  with  symmetric  power  spectra,  a  meaningful  estimate  of  time  delay 
can  be  extracted  from  this  phase  information  and  knowledge  of  the  carrier 
frequency  of  the  signal. 

The  estimator  is  evaluated  while  operating  in  white  Gaussian  noise 
which  is  in  general  correlated  between  channels,  and  curves  are  given  for 
the  density  function,  mean,  and  variance  of  the  estimator  for  various 
noise  assumptions.  The  estimator  is  shown  to  take  advantage  of  the 
processing  gain  of  large  time-bandwidth  product  signals  to  reduce  the 
variance  of  the  time  delay  estimate. 
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CHAPTER  1 
INTRODUCTION 

The  problem  of  time  delay  estimation  has  received  much  attention  in 
the  llterature^"^.  In  its  simplest  terms,  the  problem  is  to  estimate 
the  time  difference  of  arrival  of  similar  signals  in  two  different 
channels.  In  general,  the  literature  can  be  divided  into  two  distinct 
categories  passive  and  echo  location. 

In  the  passive  mode,  the  receiver  "listens"  to  a  source  in  each  of 
two  receiver  channels,  and  estimates  the  time  difference  between  the  two 
channels.  In  this  case,  very  little  may  be  known  about  the  form  of  the 
signal.  In  echo  location,  a  signal  is  transmitted  in  the  channels,  and 
the  receiver  "listens"  for  reflections  in  the  channel.  Here  the  signal 
form  is  known  to  be  a  (possibly  distorted)  time  delayed  replica  of  the 
transmitted  signal,  and  the  task  is  to  measure  the  difference  of  the 
time  delays  in  the  channels. 

In  each  of  these  modes  of  operation  there  are  two  common  methods 
used  to  estimate  the  time  delay:  cross-correlation  methods,  and  phase 
comparison  methods.  The  first  of  these  methods  performs  a  cross¬ 
correlation  of  the  two  received  signals,  selecting  as  the  estimate  of 
the  time  delay  that  value  which  maximizes  the  magnitude  of  the  cross¬ 
correlation. 

In  the  phase  comparison  method,  the  analytic  signal  from  one 
channel  is  conjugated  and  multiplied  by  the  signal  from  the  other 
channel  and  the  phase  of  the  resultant  product  is  averaged  over  the  time 
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duration  of  the  signal.  In  order  to  assign  a  meaningful  estimate  of  the 
time  delay  from  the  phase  Information,  the  signal  must  be  assumed 
narrowband. 

In  the  literature  cited  above,  these  estimators  and  variants 
thereof  have  been  extensively  analyzed  under  various  hypotheses. 

This  thesis  presents  a  method  of  time  delay  estimation  which  was 
originally  developed  by  Ricker^.  This  method  is  inherently  an  echo 
location  estimator  of  time  delay  which  uses  a  comparison  of  the  phase  of 
the  matched  filter  for  each  channel  to  estimate  the  delay.  It  will  be 
seen  that,  under  the  proper  assumptions,  this  estimator  is  capable  of 
using  large  time-bandwidth  product  signals  and  still  giving  a  meaningful 
estimate  of  the  time  delay  from  the  phase  information. 

Other  advantages  of  this  method  include  the  fact  that  it  handles 
Doppler-shifted  channels  with  ease,  and  that  with  proper  signal  design, 
it  can  isolate  multiple  scatterers  in  the  channels  and  estimate  the  time 
difference  for  each  scatterer.  With  this  estimator,  one  may  utilize 
signals  with  large  time-bandwidth  products  that  resolve  well  both  in 
time  and  in  frequency  so  that  one  can  simultaneously  estimate  the  total 
propagation  delay,  the  time  stretch  (Doppler),  and  the  time  difference 
of  arrival. 

The  next  chapter  will  present  the  receiver  structure  and  the 
estimation  procedure.  Chapters  3,  4  and  5  will  evaluate  the  perfor¬ 
mance  of  the  receiver  operating  in  additive  white  Gaussian  noise,  and 
Chapter  6  will  summarize  and  suggest  further  work  in  this  area. 
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CHAPTER  2 

RECEIVER  STRUCTURE 

This  chapter  will  discuss  the  method  in  which  the  received  signals 
are  processed,  and  will  present  notation  and  assumptions  that  will  be 
used  in  the  chapters  that  follow.  The  receiver  structure  is  shown  in 
Figure  2-1.  It  is  assumed  that  a  signal  f(t)e-^^o*',  was  transmitted 
and  that  the  received  signals,  and  r2(tj^),  are  samples  of  two 

time-delayed,  time-stretched  replicas  of  the  transmitted  signal  plus 
additive  noise.  The  function  f(t)  is  known  as  the  complex  envelope  of 
the  signal,  and  is  equal  to  2TTf^,  where  f^  is  the  carrier  frequency. 
The  received  signals,  rj^  and  r2  can  be  written  as  follows 

_ _  -  ,  s(  t .-  T,  ) 

r^(t^)  =  /e^  f[s(t^-Tj)]e  o  1  1  +  n^(t^) 

_ _  r  ,  iu>  s(t  -  T.N 

r2(ti)  =  /e^  f[s(t^-T2)]e  o  i  2  +  n2(t^) 

It  is  assumed  both  signals  are  of  the  same  energy,  and  that  the  complex 
envelope  is  normalized  such  that 

N 

1=1 

so  that  in  the  absence  of  noise,  the  energy  of  each  sampled  signal  is  E^. 
The  time-stretch  factor,  s,  is  assumed  to  be  the  same  in  each  channel, 
and  is  related  to  the  Doppler  shift,  !*>d ,  as 


f  [s(ti-Tj^)] 


k  =  1,2 


(2-1) 


(2-2) 


Gce 
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The  time  delays  for  each  channel  are  denoted  by  and  ^2*  bhe  dif¬ 
ference  in  time  delays,  At  =  -  Tp  is  the  quantity  to  be  estimated. 

The  noises  n^(t^)  and  assumed  to  be  white,  complex  proc¬ 

esses  with  circular  symmetry8»9^  snd  with  variance  of  real  and  imaginary 
2 

parts  equal  to  o^.  The  real  and  imaginary  components  for  each  channel 
are  denoted  as 

=  xj^(ti)  +  Jyj^(ti)  k  =  1,2. 

The  received  signals  are  processed  according  to  Figure  2-1.  The 
processing  signal,  r(tj^),  is  a  time-delayed,  time-stretched  replica  of 
the  transmitted  signal  with  hypothesis  time-aelay  and  time-stretch 

A  A 

values  of  x  and  s  obtained  through  some  previous  estimation  procedure. 
The  processing  signal  is  denoted  as  follows: 


f(c^)  .  . 

The  processing  siguar  is  assumed  normalized  to  unit  energy,  l.e. 


2  ? 

r(t^) 

=  I 

i=l 

f Is(t^-T)] 

is 


In  the  absence  of  noise,  the  output  of  each  summation  in  Figure  2-1 


=  /E  I  f*(s(t  -01f[s(t  -T 

O  1.  IK 


(2-3) 
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which,  when  viewed  as  a  function  of  t  and  s,  can  be  considered  as  a 
two-dimensional  correlation  between  the  received  and  processing  signals. 

The  estimates,  t  and  s,  are  often  obtained  by  processing  the 
received  signal  with  a  matched  filter.  The  matched  filter  concept  is 
based  on  the  examination  of  the  magnitude-squared  of  a  function 

A  A  A  A 

X(t,s,t,s)  over  an  appropriate  grid  of  t  and  s  values,  choosing  as  the 
estimate  of  the  time-delay  and  time-stretch  those  t,s  pairs  which 
maximize  the  magnitude  squared,  |x(t,s,t,s)|  .  Note  that  for  signals 
with  sufficient  resolution,  the  matched  fili-er  approach  is  capable  of 
resolving  multiple  reflections,  identifying  the  t  and  s  values  of  each 
reflector  as  distinct  maxima  on  the  t,s  grid. 

An  important  function  that  arises  out  of  the  discussion  of  matched- 
filters  is  the  ambiguity  function.  The  ambiguity  function  depends  only 
on  the  particular  signal  employed,  and  is  equivalent  to  the  matched  fil¬ 
ter  for  a  noiseless  scattering  channel  consisting  of  a  point  scatterer. 
The  accuracy  with  which  one  can  estimate  'Ime-delay  and  -stretch 
depends  on  the  behavior  of  the  signal  ambiguity  function  near  the 
origin.  For  the  reader  unfamiliar  with  matched  filter  concepts,  rather 
thorough  treatments  are  given  in  Van  Trees and  in  Cook,  and  Bernfeld^^ 
for  the  narrowband  case,  and  in  papers  by  Altes^2,13  for  the  wideband 
formulation . 

The  time  delay  estimator  of  Figure  2-1  works  as  follows.  First, 
the  received  signal,  r(t),  representing  a  composite  of  the  signals  r^Ct) 
and  r2(t)  (or  a  related  signal)  is  matched-f iltered  over  an  appropri- 
ate  T-s  grid,  yielding  |x(x,s,t,s)|  .  This  function  is  examined  over 
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the  grid,  with  those  peaks  that  are  above  a  pre-selected  threshold 
representing  substantial  reflections  in  the  channel  with  different  time- 

A  A 

delay  and  -stretch  values  estimated  by  the  corresponding  t  and  s  values. 

Once  a  peak  has  been  located ,  the  t  and  s  values  corresponding  to 
the  peak  are  used  to  process  the  received  signal  for  channels  1  and  2, 
yielding  Xj^(TpS,T,s)  and  X2('^2’® each  peak,  the  complex  number 
»s  ,T  ,s)  and  X2(  "^2 ^ multiplied  together,  yielding  (in  the 
absence  of  noise) 

Q  =  XjX*  =  E^*(5;f*[s(t^-T)]f [s(t^-Tp]e'^‘^o^® 

(If[s(t  -0]f*(s(t.-Tp]e-JV^-">H)ejV^^2-^l> 
i 

^  EQ*(nf*[s(tj^-T)]f  [s(t^-Tj)) 

Ij 


fls(t.-?)]f*ls(t.-T2)lej%^"-")('^i-*^j>)ejV^"2-'l^  (2-4) 

If  the  time-stretch  estimate  is  assumed  correct  (i.e  s=s)  ,  the 
first  line  of  (2-4)  reduces  to 


Q  =  EQ*(If*(s(t^-T)]f  [s(t^-Tp]  )• 


(If(s(t^-0]f*(s(t^-T2)l)ejV^ 


=  E  •(A,A-)e-^'^o 
o  12 


(2-5) 
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where 


[s(t^-T)  ]f  [s(t^-Tj^)] 


* 


=  0)  s(t  -T  ). 
0/1 


Hence,  it  is  seen  that  if  the  time-stretch  estimate  is  correct,  the  time 
difference  of  arrival  can  be  obtained  from  the  phase  of  0  as 


u  s 
o 


if  the  quantity  (Aj^A^)  is  real.  Note  that  if  =•  then  A^  =  A^  so 

•k 

that  the  quantity  Aj^A2  Is  real.  However,  with  arbitrary  At,  A^  and  A2 
must  be  real  Independent  of  one  another  for  most  signals  of  interest, 
and  it  is  desired  to  find  sufficent  conditions  for  which  this  occurs. 

To  do  this  it  is  more  convenient  to  utilize  the  continuous-time  domain. 
In  this  domain 

f  *  '•v 

Aj^  =  Jf  (st-sT)f (st-sTj^)dt, 

where  the  energy  conserving  factors  have  been  Ignored,  as  they  are 
unimportant  in  the  discussion.  Defining 

g(u)  =  f(u-a)  <=>  G((d)  =  e"^“®F(a)) 

h(u)  =  f(u-b)  <=>  HCu)  =  e~^“'’F((D), 

Parseval's  relation, 

OD  00 

r  *  r 

j  g  (u)h(u)du  =  /  G  (u))H(m)dw 


00 
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can  be  used  to  obtain 

=  /  G  (a))H((»i)dai  =  /  |F(a))|^  ^^do). 

»oo  ^co 

If  |F(-a))|  =  |F(w)|,  this  reduces  to 
00 

=  2  /  |F(a))|^  cos  t>)s(Tj^-T)da),  (2-6) 

o 

which  is  a  real  quantity.  Hence,  a  sufficient  condition  for  Aj^  to  be 

real  has  been  established.  If  the  complex  envelope  has  a  symmetric 

power  spectrum,  lF(-a))|  =  |F(a))|  ,  then  A^  is  real.  Note  that  this 

is  dependent  on  the  proper  time-stretch  estimate  but  is  Independent  of 

the  delay  estimate.  Throughout  this  thesis,  it  is  assumed  that  this 

condition  is  satisfied  so  that  Aj  and  A2  will  be  real  quantities. 

The  performance  of  this  estimator  depends  on  two  factors.  First, 

the  additive  noise  inherent  in  the  channels  will  degrade  performance  by 

adding  unwanted  terras  to  (2-5),  thereby  corrupting  the  phase,  ^  . 

o 

A  * 

Second,  the  estimates,  t  and  s  will  not  exactly  match  the  true  values. 

A 

If  the  estimate  x  is  reasonably  close  to  the  true  values  x^  and  x^, 
performance  will  not  be  adversely  affected,  since  the  only  effect  of 
this  is  that  the  factors  Aj^  and  A2  in  (2-5)  will  become  smaller  relative 
to  the  additive  noise  terms  as  can  be  seen  from  the  cosine  term  in 
(2-6).  If  the  estimate  s  is  incorrect,  however,  the  estimator  becomes 
biased  due  to  the  additional  exponential  that  appears  in  (2-4).  The  x 
and  s  estimates  obtained  from  the  matched  filter  will  obviously  become 
less  accurate  as  the  noise  level  increases. 
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In  Chapters  3  and  4,  the  first  of  these  problems  is  addressed,  i.e. 
that  of  the  unwanted  noise  terms.  In  these  chapters  it  is  assumed  that 
the  estimates  of  t  and  s  are  correct.  In  Chapter  5,  the  effect  of  a 
T-s  mismatch  is  discussed  using  a  Cramer-Rao  lower  bound  approach.  To 
simplify  the  analysis  and  notation,  it  is  assumed  throughout  the  follow¬ 
ing  that  the  channels  contain  only  one  scatterer. 
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CHAPTER  3 

UNCORRELATED  NOISE 

The  first  case  to  be  considered  is  that  of  the  receiver  operating 

in  white  Gaussian  noise  uncorrelated  between  channels.  The  complex 

noise  process  is  assumed  stationary,  ergodic,  zero  mean,  circularly 

2 

symmetric,  with  variance  of  real  and  imaginary  parts  equal  to  o^.  The 
real  and  Imaginary  components  for  each  channel  are  independent ,  and 
are  denoted  as  follows: 

niCtf)  =  Xj(t^)  + 


for  all  i,j.  (3-1) 
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The  magnitude-squared  of  the  noise  process  has  expectation 

E(|n^j2)  .  *  y2.)  -  2,2  k.1.2. 

The  received  signal  for  each  channel  at  time  t^  is 


or 


^  f[s(t.-Tj^)]ej“o®^4  V  +  n^(t.) 


r.  .  =  /^  f.  .  ej“os(t^-T^)  ^ 
ki  o  ki  Tci 


where  f^^^  i  f[s(t^-Tj^)]  is  the  complex  envelope  of  the  received  signal. 
If  the  processing  signal  f ( t .  )  is  denoted  as  follows: 


Ht.)  =  g-ja)oS(ti-0 


where  f^  =  f  [s(t^-T)]  represents  the  complex  envelope  of  the  processing 
signal,  then  the  output  of  each  summation  block  in  Figure  2-1  can  be 
written 


v=rr,-  -rf  =“j“^s(t.-T) 

\  j:  i  ki  ^  i  ^ 


/e;  VV  . 


/g-  f  f  ej‘^oKs-s)t^+ST-sT,^]  f  -ji^osCti-T) 

o  i  ki  ^i  i 


/T'f.f,.  +  It  = 

o  1  ki  ki  1 


(if  s=s)  (3-2) 


where  the  last  line  assumes  the  proper  time-stretch  estimate,  s  =  s. 
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Since  the  Gaussian  noise  components  are  all  assumed  independent  of 
one  another,  the  explicit  and  implicit  phase  factors  in  the  second  term 
of  (3-2)  can  be  absorbed  into  the  noise  process  without  altering  the 
statistics  of  the  situation.  This  yields 

Xfc  ■  .  I  k  .  1.2.  (3-3) 

The  second  term,  being  the  sura  of  independent,  zero-mean  Gaussians  is 
itself  a  zero-mean  Gaussian,  which  shall  be  denoted  by  The  real  anc 

Imaginary  parts  of  this  term  have  variance 


.2  2 


=  I  |f  1^  =«  0^  , 

n  J  '  i '  n  ’ 


since  Z  |fj^|  =  1 ,  so  that  (3-3)  can  be  written  as 


(j^  =  \  k  =  1,2, 


where  =  E  f^fj^^. 

The  output,  Q,  of  the  estimator  of  Figure  2-1  is  then 


*  r;r-  ,  Im  s(t-T,)  , -  -lu)  s(t-T„)  * 

Q  =  XiXo  “  o  1  +  n,  •  /E  A^e  o  2  + 

12oi  loz  Z 


A,A.,  e^V^^2“''l^  +  /e"  A,n.  +  /T  A. n!  +  n,n* 


o  1"2 


o  2  ’1  o  12  '1  7 


(3-4) 


where,  as  above,  the  complex  exponentials  have  been  absorbed  into  the 
noise  terms  without  altering  the  statistics. 

It  has  been  previously  assumed  that  s=s.  If  it  is  now  also  assumed 


that  the  time  difference  of  arrival.  At,  Is  small  and  that  the  time  delay 
estimate,  t  Is  accurate  such  that  t  "  Tj^  (k  =  1,2),  then  fj^^  •  f^^, 
so  that  Aj^  •»  Z  |?j^|  =  1. 
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To  see  what  Is  meant  by  the  words  "small"  and  "accurate,"  one  needs  to 
examine  (2-6)  where  it  is  seen  that  if  cos  o)  s(t,-t)  “>  1,  then  A.  will 

1C  iC 

be  at  its  maximum,  which  is  1.  Here  to  denotes  the  maximum  frequency 

max  ^  ^ 

of  the  complex  envelope  (i.e.,  equals  one-half  the  signal  band¬ 

width).  In  this  case  (3-4)  reduces  to 


Q  =  E  e 
o 


jaioS(T2-Ti)  ^  ^  ^  ^  ^  ^  * 

o  1  o  2  12 


(3-5) 


The  presence  of  the  last  term  in  (3-5)  renders  further  analysis 
intractable,  so  it  is  desired  to  quantify  conditions  under  which  this 
term  becomes  negligible.  The  expectation  of  the  magnitude-squared  of 
of  the  process  can  be  found  as 

Edn^n*!^}  =  E{[Re(n^n*)]^}  +  E{[im(n^n*)l^} 

=  4a  , 
n 

where  E{*}  is  the  expectation  operator.  Thus,  if  the  input  signal-to- 
noise  ratio  is  defined  as 


h  = 


Input  Signal  Energy 


Expected  input  Noise  Energy  N.E{|n,,,|^)  2No^’ 

'  k.t '  n 


then  the  ratio  of  the  magnitude-squared  of  the  first  term  of  (3-5)  to 
that  of  the  second  or  third  term  has  expectation 

E^  E 

- - - IT-  =  -4-  =  Nh  (3-6) 


E{Eo*l\l  > 


2o 


while  the  ratio  of  the  magnitude  of  the  first  term  to  that  of  the  last 


term  has  expectation 
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=  (Nh)^  . 
4  a 


(3-7) 


Here,  as  In  the  previous  chapter,  N  Is  the  number  of  input  samples.  By 
comparing  (3-6)  and  (3-7),  it  is  seen  that  the  last  term  of  (3-5)  is  of 
second  order.  For  example,  if  there  is  a  OdB  input  SNR  (h=l)  and  there 
are  1000  input  samples,  then  the  second  and  third  terras  are  30dB  below 
the  first,  while  the  last  term  is  60dB  below  the  first.  For  the  rest  of 
the  development,  it  is  assumed  that  the  quantity  Nh  is  large  enough  so 
that  the  last  term  of  (3-5)  can  be  Ignored.  It  will  be  seen  that  for 
values  of  Nh  as  small  as  10,  this  approximation  yields  variances  con¬ 
sistent  with  those  obtained  through  computer  simulations. 

It  is  seen  that  for  large  Nh,  the  output  SNR,  given  approximately 
by  (3-6),  is  N  times  as  large  as  the  input  SNR.  The  factor  N  can  thus 
be  considered  as  a  processing  gain.  For  a  receiver  sampling  at  a 
frequency  equal  to  the  bandwidth  of  the  signal,  N  is  equal  to  the 
time-bandwidth  product  of  the  signal,  so  that  the  receiver  structure  of 
Figure  2-1  is  seen  to  have  a  processing  gain  equal  to  the  time-bandwidth 
product  of  the  signal. 

Defining  such  that 

_  * 

(nj+n2), 


(3-5)  reduces  to 


Q  =  V'P  +  n^. 


(3-8) 


where  is  a  complex  Gaussian  process  with  variance  of  real  and 
imaginary  parts  equal  to 


2  oc  ^2 
o.  =  2E  a 

3  on 


(3-9) 
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Q  can  be  written  in  terms  of  its  real  and  imaginary  parts  as 

Q  =  E 

=  E  (cosd)  +  isind)  )  +  u+iv 
=  a  +  u  +  j(b+v) 


where 


V''2-b> 


a  = 


E  cosdi  ,  b  =  E  sind) 
o  o  o  o 


u  = 


Re{n2},  var(u)  = 


V  = 


Im{ri2},  var(v)  =  , 


so  that  our  estimate  of  denoted  by  (J) ,  is 


In  the  absence  of  noise  it  is  noticed  that  as  desired. 

In  order  to  evaluate  the  performance  of  this  estimator,  it  is 
desired  to  find  the  probability  density  function  for  <j) ,  denoted  by  f(4>)» 
The  density  function,  f(<j))  would  be  more  precisely  notated  as  f(4)|4i^), 
and  might  be  more  properly  regarded  as  a  likelihood  function.  Also, 
it  is  noted  that  the  estimator  to  be  derived  is  a  maximum  likelihood 
estimator^  1 » ^5  of  for  the  uncorrelated  noise  case. 

The  numerator  and  denominator  of  the  arctan  argument  in  (3-10) 
are  Independent  and  have  Gaussian  densities  with  mean  b  and  a, 
respectively.  The  joint  density  for  numerator  and  denominator  is  thus 
the  product  of  the  two  marginal  Gaussian  densities,  and  to  find  the 
density  for  (ft,  one  may  convert  to  polar  coordinates  as  follows,  and 
Integrate  over  0  <  r<  ®.  Letting  u+a  =  rcosij)  and  v+b  =  rslnij),  the  joint 


density  becomes 
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f(u,v)  =  — exp  — ^  (u^+v^) 


2x0: 


2o: 


f(r,4>)  =  — r  exp  — ^ 


2  2 
[(rcos<J>-a)  +  (rsin!j)-b)  ]  , 


2xa: 


2a, 


where  the  factor  r  is  the  Jacobian  of  the  transformation.  The  desired 
density  is  obtained  by  integrating  over  r: 

CO 

f(<{>)  =•  /  f(r,4')dr 

o 


exp 


lira: 


2a: 


[a  +b  ]  *  /  r  exp  — ^  [r  -  2r(acosij)  +  bsin(j»)Idr. 


2o: 


Recalling  that  a  =  E  cosij)  and  b  =  E  sini}>  ,  and  using  a  well  known 

o  o  o  o’ 

trigonometric  identify,  the  above  expression  becomes 

,2  „ 


f(^) 


1 


-E 


2xa: 


exp 


20^  o 


/  r  exp  — ”  2E^^cos(<>-|^>^)]dr , 
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or,  making  a  change  of  integration  variable, 

2  2 
E  <»  "^9 

f(4i)  =  - 1—  /  r  exp  — ^  [r  -  2  cos(4)-iji^)r  +  l]dr. 


2x0^  o 


2a: 


which  is  a  function  only  of  <(>~t 


h  =  E  /2Na  so  that 
o  n 


2  2  2 
and  f — )  .  Recall  that  o  =  2E  o  ,  and 
a.  3  o  n 


2E^a 
o  n 


=  Nh. 


Here  again,  h  is  the  input  signal  to  noise  ratio,  and  N  is  the  number  of 
input  samples.  Finally,  then,  the  density  is  written 
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f('l')  =  e  2  /  r  exp  — ^  -  2  cos((j)-(ti  )r]dr,  (3-11) 


-Nh  “> 


-Nh  ,  2 


2it 


which  is  a  function  only  of  snd  the  quantity  Nh. 

Unfortunately,  the  density  given  by  (3-11)  cannot  be  expressed  in 
closed  form.  The  integral  in  (3-11)  can  be  found  in  the  tables  by 
Gradshteyn  and  Ryzhik^^  (//3. 462-5),  where  it  is  expressed  in  terms  of 
the  error  function,  but  the  solution  is  only  valid  in  our  case  for 
cos( i|i-(j)^)  <  0.  The  density  of  (3-11)  is  Integrated  nuiiicrically  on  a 
VAX  11/782,  yielding  the  curves  in  Figures  3-1  and  3-2.  Note  that 
although  these  densities  are  defined  over  -it  <  it'-ii'o  ^  tt  ,  in  Figures  3-1 
and  3-2  they  are  shown  over  a  smaller  range  to  better  see  their  shape. 
Also,  a  curve  of  the  variance  as  a  function  of  the  output  signal  to 
noise  ratio  Nh  was  generated  digitally,  yielding  Figure  3-3. 

It  must  be  noted  that  the  moment  quantities  depend  on  the  inte¬ 
gration  limits  chosen.  Mathematically,  any  limits  of  the  form 
0  <  (fi  <  9  +  2Tr  will  do,  however  for  the  application  here,  the  choice 
is  clear.  Choosing  the  limits  <}'q-tt  <  4>  <  ‘J’o'*'''' «  mean,  y,  is  given  by 


dl  +TT 

o 


TT 


u  =  /  ii>f(i>)di})  =  /  (<ji+({>^)f (di+ij>^)d'j) 

(p  -TT  -IT 

O 


IT  tr 

=  /  <{)f ( i)+<|>^)d*>  +  /  f ( i{>+iJiQ)d4> 

-IT  -TT 


=  0  +  1^  =  1 
o  o 


Nh=100 
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IKUOQ 


Figure  3-1:  Probability  Density  of  Estimator  for  Uncorrelated  Noise;  Nh=10, 20, 50, 100 


Probability  Density  of  Estimator  for  Uncorrelated  Noise;  Nh=100, 200, 500, 1000 
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(pdjonbs  s0Oj6»p)  ooudijoa 


no^ 
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where  the  first  Integral  is  zero  due  to  its  odd  integrand  integrated 
over  symmetric  limits,  and  the  second  integral  is  1  since  it  represents 
the  area  under  density,  fC’t')*  Hence,  with  this  choice  of  integration 
limits,  the  mean  is  equal  to  the  true  value  (i.e. ,  the  estimator  is 
unbiased) ,  so  that  the  variance  is  a  measure  of  the  fluctuation  of  i|i 
about  the  true  value, 

The  estimate  of  time  delay  is  given  by 


(0  s 
o 


which  has  mean  equal  to  the  true  value  (if  s=s),  and  variance 


E{At}  = 


u 


(0  s 
o 


bi  S 

o 


var{AT}  =  — i — ~  , 

o 


where  o  =  var(<{))  is  shown  in  Figure  3-3. 


From  Figure  3-3,  it  is  seen  that  the  variance  of  this  estimator 
decreases  approximately  as  1/Nh.  Under  previously  stated  assumption,  N 
can  be  replaced  by  the  time-bandwidth  product  of  the  signal,  and  the 
time  delay  estimator  is  seen  to  take  advantage  of  the  processing  gain  of 
large  time-bandwidth  product  signals  to  reduce  the  variance  of  the  time 


difference  estimate. 
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In  order  to  test  the  validity  of  the  variance  curve  of  Figure  3-3, 
the  estimation  procedure  of  Figure  2-1  is  simulated  on  a  computer.  The 
signal  used  for  the  simulation  is  based  on  a  Welsh  construction^^  with  a 


real  envelope.  The  signal  is  of  the  following  form 


f(t)ej 


2TTf  t 
o 


where 

f  =  30,000Hz 
o 

f(t)  =  w(  t)cos2Tif^t  (i-l)Tg  <  t  <  ITg 

Tg  =  subpulse  length 
=  .05  seconds 

w(t)  =  50  dB  Taylor  Window^^ 


i 
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11 
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12 
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The  ambiguity  function  of  this  signal  is  shown  in  Figure  3-4.  Note 
that  this  signal  has  good  resolution  both  in  t  and  in  s,  due  to  its 
large  time- band width  product.  Note  also  that  by  giving  the  signal  a 
real  envelope,  one  obtains  a  symmetric  spectrum  as  required  by  (2-6). 

The  signal  is  sampled  at  1700Hz  and  matched-f iltered  using  the 
correct  t  amd  s  values,  t  =  t ^  =  1 .0,  s  =  s  =  1 .  Next,  1000  complex 
Gaussian  noise  samples  were  generated,  normalized  as  in  (3-9),  and  added 


COSTAS  CODE:  TAYLOR  WINDOW 

REAL  ENVELOPE  ON  TRANSMITTED  SIGNAL 


Figure  3-4:  Ambiguity  Function  of  Signal  Used  in  Example 
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to  0  as  in  (3-8)  to  obtain  a  sample  mean  and  variance  to  compare  with 

the  theoretical  unbiased  mean  and  the  variance  of  Figure  3-3.  The 

results  are  summarized  in  Table  3-1  table  for  =0  and  <Ji  =  45°. 

o  o 

Table  3-1 

Theoretical  versus  Sample  Means  and  Variances 


Nh 

10 

100 

1000 

10,000 

1 

iTheoretical  Mean  (Deg) 

0.000 

0.000 

0.000 

0.000 

=0°  [Sample  Mean 
o 

0.068 

0.045 

0.015 

0.005 

[Theoretical  Variance 

1 

380.0 

33.2 

3.3 

0.36 

1 

[Sample  Variance 

1 

369.0 

34.0 

3.4 

0.34 

[Theoretical  Mean 

1 

45.00 

45.00 

45.00 

45.00 

1 

=  45°  [Sample  Mean 
o  1 

45.84 

45.25 

45.09 

45.05 

1 

[Theoretical  Variance 

1 

380.0 

33.2 

3.3 

0.36 

I 

[Sample  Variance 

379.6 

1  35.0 

3.5 

0.34 

It  is  seen  that  the  expressions  derived  above  for  the  density  func¬ 
tion  of  the  estimator  generally  predict  variances  consistent  with  the 
sample  variances  of  the  simulations  even  for  values  of  Nh  as  low  as  10, 
and  the  sample  mean  approaches  the  theoretical  mean,  i}>^,  for  large  values 
of  Nh. 

This  concludes  the  analysis  of  the  estimator  operating  in  uncorre¬ 
lated  noise.  The  next  chapter  will  discuss  the  considerably  more  com¬ 
plicated  case  of  noise  which  is  correlated  between  channels. 
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CHAPTER  4 
CORRELATED  NOISE 


In  this  chapter,  the  receiver  performance  is  evaluated  when 
operating  in  white  Gaussian  noise  which  is  correlated  between  channels. 
As  in  the  previous  chapter,  the  noise  is  assumed  to  be  a  white,  cir¬ 
cularly  complex  Gaussian  process,  so  that 


^^’^ki^ki^  =  0  for  all  i;  k  =  1,2 


=  ^tyki^kj^  =  %  ^ij  =  ^*2. 


The  noise  process  is  assumed  to  be  correlated  as  follows:^ 


-  E(y,iy2j)  - 


E{Xjjy2j)  -  -Efxj^yn)  =  A<,_,  . 


A  physical  interpretation  of  p  and  X  is  given  in  Appendix  A.  It  is 
noted  here  for  future  reference  that  p  and  X  are  such  that 


9 

p  '  +  X 


1. 


(4-1) 


This  can  be  seen  from  examining  eigenvalues  of  the  covariance  matrix 

for  the  density  f (XpypX2  ,y2) »  noting  that  since  the  covariance  matrix 
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is  positive  definite,  its  eigenvalues  must  all  be  positive. 

The  development  of  Chapter  3  may  be  followed  up  to  (3-2)  without 
modification.  Equation  (3-2)  is  repeated  here  for  convenience: 


k  =  1,2. 

(4-2) 
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represents  the  output  of  the  summation  block  in  Figure  2-1 ,  and 
(4-2)  assumes  the  correct  tine-stretch  estimate.  Since  the  complex 
noise  is  assumed  correlated  between  channels,  more  care  must  be  taken 
when  absorbing  the  complex  exponential  of  the  last  terra  of  (4-2)  into 
the  noise  process. 

If  the  second  term  of  (4-2)  is  written  in  terms  of  its  real  and 
imaginary  components, 


E 

i 


-jtil  s(t 
e  o 


“k 


then  the  correlation  coefficients  for  Uj ,  ,  U2 ,  can  be  expressed 

in  terms  of  the  input  correlation  coefficients,  p  and  X.  Using  the 
results  of  Appendix  B  and  the  fact  that  the  time  samples  of  the  noise 
are  independent,  it  can  be  shown  that  ,  Vj ,  U2,  V2  have  the 
same  correlation  coefficients  as  ’^21*  ^2i’ 

2 

E{u,u-}  =  E{v  V  }  =  po 
12  12  n 

2 

E{u,v„}  =  -E{u_v  }  =  Xa 
12  2  1  n 


It  is  noted  that  the  variances  also  remain  unchanged,  due  tc  the  fact 
that  the  processing  signal  is  of  unit  energy. 

As  in  the  previous  chapter,  this  processed  noise  term  is  denoted 


by 

"k  ■  “k 

with  variance  of  real  and 


■  ?  "ki  h 


imaginary  parts 


A  A 

S(t^-T) 

equal  to  o 


2 

n 


k  =  1,2 
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If  is  again  defined  as 


I 


k  =  1,2, 


then  the  output  of  each  summation  block  in  Figure  2-1  can  be  written  as 


=  /e~  a,  e^“o®^^~^k^  +  n,  k  =  1,2 
ok  k  ’ 


so  that  the  output  Q,  is  given  by 


Q  -  X^X*  =  [/TAj  +  n*l 


-  E  •a,a,a^V<^2‘'i'  +  /rA,n,e'jV<"‘’2'  +  /Ta, n*A^ V"‘^l ’ 
ol2  o21  ol2 


+  n,n* 


(4-3) 


Defining 


"l  = 


„  -ju  s(t-t  )  -  - 

n^e  J  o  2  =  u^  +  jv^ 


^2 


„  -io)  s(t-t  )  -  - 

n2e  o  1  =  u^  +  jv^ 


it  is  again  desired  to  find  the  correlation  coefficients  of  u^^ , 

U2 ,  '^2  terms  of  P  and  A,  The  results  follow  from  direct  application 
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of  Appendix  B,  and  are  summarized  below: 


=  ECv^v^}  =  (cp  -  dX)o 


2 

n 


where 


E{u  V  }  =  -E{v  u  }  =  (cX  +  dp)a^  =  Xa^ 
12  12  n  n 

A 

c  =  cos  m  s( T.-T, ) 
o  2  1 

d  =  sin  u)^s(T2-Tp. 


and  the  new  correlation  coefficients,  p  and  X  are  introduced  to  simplify 
notation  in  the  following  development.  It  is  noted  for  future  reference 
that 


;;2  .  r2  _  (1  ,  .2w^2  .  .1. 
p  +  X  =  (c  +  d  )(p  +  X  ) 


2  2 

=  p^  +  X^  <  1 


where  p^  +  X^  <  1  as  in  (4-1). 


Equation  (4-3)  can  now  be  written  in  terms  of  and  n2  as 


Q  =  E  .A,A_e^V^^2'^^  +  fTk'rx.  +  /Ta,  n*  +  n,  V^^2‘’'l^ 
012  o21  ol2  12 


=  E^*e^‘^o®^^2  ^1^  +  "^(n^  +  02  e^“o^^^2  ^1^  ,  (4-4) 


where  the  last  step  assumes  that  the  estimate  of  time  delay  is  correct, 
and  that  the  time  difference  of  arrival  is  small  so  that  A^  “  1  and 
A2  “  1 ,  as  in  the  previous  chapter. 
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Again,  as  a  concession  to  tractability ,  conditions  under  which  the 
last  term  is  negligible  are  quantified.  The  expectation  of  the 
magnitude-squared  of  the  process  njn2  can  be  found  directly  through 
tedious  integration  to  be 


E{|n^n*j2}  =  E{[Re(f;^n*)]2}  +  E{[  im(n^n*)  ] 

4(1  +  <  8a^ 

n 

-2  **-2 

Where  the  inequality  arises  from  the  fact  that  p  +  X  <1.  Thus 
with  the  input  signal-to-noise  ratio  defined  as 


E 

o 

2 

2Na 

n 


> 


the  ratio  of  the  magnitude-squared  cf  the  first  terra  of  (4-4)  to  that  of 
the  second  or  third  term  has  expectation 


E((/rsj2} 


Nh 


(4-5) 


as  in  (3-6),  while  the  ratio  of  the  magnitude-squared  of  the  first  term 
to  the  expected  magnitude-squared  of  the  last  term  is 


Ertln^n^l 


■}  8a 


n 


(4-6) 


From  (4-5)  and  (4-6)  it  is  seen  that  the  last  term  of  (4-4)  is  again  of 
second  order,  becoming  negligible  for  large  values  of  Nh.  For  example, 
if  Nh  =  1000,  then  the  second  and  third  terms  are  30dB  below  the  first. 
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while  the  last  term  Is  at  least  57dB  below  the  first.  Throughout  the 
rest  of  the  chapter,  the  last  term  of  (4-4)  is  dropped. 

In  terms  of  its  real  and  Imaginary  parts,  0  can  be  written  as 


0  =  E 


•e^v^ 

)  o  1  o  2 


where 


=  E  -e^^  +  +  /^* 

o  o  1  o  2 


=  a  +  jb  +  +  U2  +  j(v^  -  v^)] 

=  a  +  /^(u,  +  u„)  +  i[b  +  /e~(v,  -  V-)] 
o  1  z  01  z 


♦o  ■  V'"2  - 


a  =  E  cos<|) 
o  o 


b  =  E  sini}) 
o  o 


u^  =  Re{nj^} 


vi  =  Im{h^} 


U2  =  Re{n2) 


=  Im{n2} 


var{u^}  = 


var{v,  }  =  a 
k  n 


-  _  -  -  -  2 

E{u,u„}  =  E{v  V  }  =  p  a 
12  12  n 


E{u,v  }  =  -E{v  u„}  =  X  o 
12  12  n 


P  =  P  cos<Ji  -  \  sinij) 
o  o 


k  =  1,2 


=  p  sinif  +  X  cosi^ 
o  o 


X 
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The  estimate  of  4*^  is  then  given  by  (^ ,  where 


=  tan  [ - 

a  +  /e~ (u  +u  ) 
o  1  2 


^  ^a  +  x'^ 


and  where  x  h  Uj^+  ,  y  =  v^-  la  order  to  evaluate  the 

statistics  of  this  estimator,  it  is  necessary  to  find  the  covariance 
matrix  for  the  numerator  and  denominator  of  the  argument  of  the  arc¬ 
tangent.  Denoting  this  by  it  is  seen  that 


R 

-xy 


E{x^}  E{xy} 

2 

=  2E  -0 

(1  +  P)  -  X 

^{xy}  E{y^2 

o  n 

-X  (1  -  Pi 

o  poo 

X  xy  x  y 

poo  o  2 
xy  X  y  y 


where 


=  2E  +  p) 

X  on 


=  2E  -  p) 

y  on 


-X 


xy 


/l-p2 


-2  ~2 

It  is  noted  that  p  <  I  since  p  +  X  <1.  The  joint  Gaussian  density 
xy 


is  then 


8 


f(x.y)  =  - - rTy  \  ^ 

2r|R  2 

I 


“-xy 


X 

y 
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1  -1 
- — -  exp  - r~S-T' 

ZifO  a  /l-p^  °x 

X  y  xy  —  ^ 


-^»xy  ^  ^  VI 

X  y  a 

y 


Once  again,  the  change  of  variables  x  +  a  =  r  cosi|i,  y  +  b  =  r  sin^i  is 
made  to  obtain  the  density  for  <p  =  tan  — J.  After  considerable 

algebraic  manipulation  the  joint  density  f(r,i}i)  is  obtained: 


where 


f(r,(j))  = 


-k  +  2v  (— ) 

e  e  '■E  '■E 

o  o 


2ro  a  /l-p 
X  y  xy 


2  2 

-  1  3  ^  ^  o  ab 

2  'I  2  2  ^xy  o  a 


E  2  2 

o  cos  <t>  sin  <1 

2(l-p2  ) 

xy  X  y 


sin^icos^i 

D  "  ' 

xy  0  0 
X  y 


^o  a  cos(}>  .  b  sin(}) 

2  2^2 

2(l-p  )  a  o 

xy  X  y 


+  (■ 
X  V  '■ 


a  siniji  +  b  cosf^i- 


The  desired  density,  fCt}’)  is  again  obtained  by  integrating  the  joint 
density,  ftr,(}i)  over  all  values  of  r,  i.e. 


f(  t')  =  /  f(r,^))dr , 


This  integral  is  of  the  sane  form  as  that  of  (3-11),  but  with  consid¬ 
erably  more  complicated  coefficients  in  the  exponential.  These  coef¬ 


ficients  can  be  substantially  simplified  by  writing  o  ,  o  and  p  in 

/  &  X  ’  y  ^xy 
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terms  oi  cr  ,  p,  and  X  and  by  recalling  that  a  =  E  cosct  and  b  =  E  sin<t  , 
n  o  o  o  o 

yielding 

E  _ 

k  =  — n — — o  [l  -  P  cos(2(j)  )  +  2  X  cos<))  sin<J>  ] 

,  2,,”2r2v^  o  o  o 

4a  (1-p  -X  ) 
n 


E  ^  ^ 

V  =  — 5 — ^-5 — j  [1  -  P  cos(2<(>)  +  2  X  cos^i  sin<Ji] 

4a^(l-p  -X^) 
n 


E 

— T — — T  [cos(<{'“'l'  )  -  P  cos(())+(|>  )  d  X  sin((J)+<J)  )].  (4-7) 
,  2 ,,  “Z  tIs  o  o  o  ^ 

4a  (1-p  -X  ) 
n 


The  density  is  then  written 

f(^)  =  - ® -  ■  /  r  exp[-u(|-)^  +  2v(I-)]dr 

,  ^  2  ,,  -2  r2  o  o  o 

4TrE  a  /1-p  -X 
o  n 

-k 

E  e 

=  - ^  —  J  j-  exp[-u  r^  +  2vr]dr  (4-8) 

/  2  ,,  ”2  -2  o 

4Tra  /1-p  -X 
n 


Where  the  second  line  results  from  a  change  of  integration  variable. 

Once  again  the  integral  for  f(<{>)  cannot  be  put  in  closed  form.  It  is 

seen  that  for  the  correlated  noise  case,  the  density  f(<{))  is  a  rather 

2  -  - 

unwieldy  function  of  several  parameters:  E  /a  ,  p  X  and  (ji  .  In 

on  p 

addition,  P  and  X  are  in  turn  functions  of  and  the  input  noise 

2 

correlation  coefficients  P  and  X.  The  Quantity  E  /a  in  (4-7)  and  (4-8) 

o  n 
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is  replaced  by  the  quantity  2Nh  to  get  the  density  in  terms  of  the  input 
signal-to-noise  ratio. 

2 

Roughly  speaking,  the  quantity  controls  the  sharpness  of  the 

density,  while  the  correlation  coefficients  alter  its  shape.  The 
density  is  integrated  numerically  for  various  combinations  of  parameters 
yielding  Figures  A-1 ,  A,  7,  and  10. 

The  densities  are  then  used  to  generate  curves  of  the  mean  and 

variance  of  the  estimator  for  these  combinations  of  parameters,  yielding 

the  continuous  curves  in  Figures  A-2,  5,  8,  11  and  A-3 ,  6,  9,  12.  It  is 

again  noted  that  the  moment  quantities  depend  on  the  Integration  limits 

chosen.  For  the  correlated  noise  case,  the  choice  of  these  limits  does 

not  seem  as  clear  as  it  was  for  the  uncorrelated  noise  case,  because  the 

estimator  is  now  biased.  If  the  same  limits  are  chosen  for  this  case 

however,  the  estimator  will  become  unbiased  for  large  values  of  Nh,  so 

the  integration  limits  are  again  chosen  as<>  -  +ir. 

o  o 

Tnese  curves  are  again  verified  using  the  same  signal  and  sampling 
parameters  as  in  the  previous  chapter,  yielding  the  discrete  points  in 
the  mean  and  variance  curves.  The  noise  model  used  in  these  simulations 
is  described  in  Appendix  A,  with  the  incoherent  noise  correlation 
coefficient,  set  equal  to  zero  so  that  the  incoherent  noise  is  uncor¬ 

related  between  channels.  The  coherent  noise  is  the  same  in  each  channel 
except  for  a  phase  delay. 


c  n 


c 

li 
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A) [ su^a 


Figure  4-1:  Probability  Density  of  Estimator; 
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Figure  A-2:  Mean  of  Estimator;  p=.67,  X-.67 
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ce 


Nh=100 
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Figure  4-4:  Probability  Density  of  Estimator;  p=0 


(p»jonbi  »9uopDA 


Nh=100 
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Figure  4-10:  Probability  Density  of  Estimator; 
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ean  of  Est 
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(p»jonbi  »9uot4DA 


k. 


Figure  4-12:  Variance  of  Estimator; 
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where  <t>  =  tan  P  and  X 

c  p 

cents.  The  incoherent  noise 


being  the  input  noise  correlation  coeffi- 
has  variance  of  real  and  imaginary  parts 


while  the  coherent  noise  has  variance 


> 


where  r  =  /p  +X  and  o  is  the  total  noise  variance. 

’  n 

From  these  figures,  it  is  seen  that  both  the  sample  and  theoretical 
means  approach  the  true  values  for  large  values  of  Nh.  Also,  the  vari¬ 
ances  predicted  by  theory  are  consistent  with  the  sample  variances  gen¬ 
erated  by  the  simulation  for  large  Nh.  The  theoretical  curves  are 
accurate  only  for  Nh  of  about  100  or  greater  due  to  the  approximation 
made  by  developing  the  last  term  of  4-4, 

Again  it  is  seen  that  for  large  Nh,  the  variance  drops  of  as  1/Nh. 
Replacing  N  with  the  time-bandwidth  product  of  the  signal  (as  in  the 
previous  chapter),  it  is  again  seen  that  the  estimator  takes  advantage 
of  the  processing  gain  of  large  time-bandwidth  signals  to  reduce  the 
estimator  variance. 

The  analysis  of  this  and  the  previous  chapters  has  assumed  that  the 
time-stretch  and  -delay  estimates  are  correct.  In  most  cases  these 
estimates  will  not  be  exactly  correct.  The  effect  of  a  x-s  mismatch  is 
the  subject  of  the  following  chapter. 
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CHAPTER  5 


T-s  MISMATCH 


In  Chapters  3  and  4  it  was  assumed  that  the  estimates  of  time-delay 
and  -stretch  were  correct.  In  practice,  these  estimates  will  not  be 
exact,  and  this  will  affect  the  performance  of  the  time  difference 
estimator.  In  order  to  examine  the  effects  of  a  x-s  mismatch,  it  is 
necessary  to  go  back  to  Equation  (3-2),  before  any  assumptions  about 

A  A 

T  and  s  were  made.  The  second  line  of  (3-2)  is  repeated  here  for 
convenience : 


-  /r.  r  f.f,,  +  £  n,  .  t 


'k  o  “  i  ki 


ki  "i 


where  represents  the  output  of  the  k  channel  in  Figure  2-1.  This 
can  be  rewritten  in  the  following  form: 


A  A  A  A 

/=—  ju)  s(t-T,  )  -jo)  (s-s)T  r  c  c  (s-s)t,  .  „ 

=  /E  e-^  o  k  e  o  •  I  f^f,  ^  e-*  o  i  +  n, 


"k  o 


=  /E~e^'^k  e^^  E  f^f.  .  e^'^i  + 


where 


♦k  ' 


O.  =  (D  (s-s)t.  ,  6  =  w  (s-s)t 

i  o  1  o 


_  ?  -jw  s(t,-T) 

\  J  "ki  'i  '  °  ‘ 
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(5-3) 
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so  that  there  is  no  bias  in  the  estimate  of  (Ji  due  to  an  incorrect  time- 

o 

stretch  hypothesis.  There  still  remains,  however,  a  bias  in  the  time 
difference  estimate.  This  bias  occurs  from  using  the  estimated  time 

A 

stretch,  s,  rather  than  the  true  time-stretch,  s.  The  estimated  time 
difference  (in  the  absence  of  noise)  is 

- (1  -  izi)  =  At(1  -  (5-4) 

m  a  m  s  +  0)  (s-a)  V 

o  o  o 

where  At  is  the  estimate  of  true  value,  4)^  =  “ 

a)^s(T2-T^) ,  and  where  a  binomial  expansion  yields  the  approximation. 

A 

For  reasonably  good  time-stretch  estimates  (i.e.,  |s-s|<<l)  this  bias 
will  be  negligible. 

Comparing  the  magnitude-squared  factor  in  (5-3)  with  (2-3)  and 
noting  the  discussion  following  (2-3)  it  is  seen  that  this  factor  is  the 
magnitude-squared  of  the  matched  filter  output  whose  main  peak,  is  at 

A  A  A  A 

T  =  Tj  ,  s  ==  s,  evaluated  at  the  point  (t,s)  corresponding  to  the  esti- 
mates  of  the  the  true  values.  If  the  T  and  s  estimates  are  not  exactly 
correct,  the  magnitude-squared  factor  in  (5-3)  will  become  smaller, 
effectively  lowering  the  signal-to-noise  ratio. 

Thus,  it  is  seen  that  a  t-s  mismatch  does  not  affect  the  noise  terms 
(or  change  any  bias  due  to  the  noise  terms),  and  if  the  time  difference 
of  arrival  is  small,  the  only  effect  is  a  slight  bias  represented  by 
(5-4),  and  a  decrease  in  the  signal-to-noise  ratio  due  to  the  magnitude- 
squared  factor  in  (5-3). 
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If  the  time  difference  of  arrival  cannot  be  assumed  small,  as  was 
necessary  to  derive  (5-3),  then  the  sum  must  remain  as  in  the  first  term 
of  (5-2).  The  problem  thus  becomes  more  difficult,  and  the  performance 
depends  on  the  particular  signal  employed. 

A  special  case  of  interest  here  is  that  of  a  real  envelope  signal 
(such  as  that  employed  in  the  examples  of  the  previous  two  chapters).  In 
this  case,  the  double  sum  in  first  term  of  (5-2)  can  be  written  as 
follows : 


I  " 


* 

I  f.f„. 

j  j 


-  £ 


.  Vii 
1  3 


^  (5-5) 


(where  the  complex  conjugate  symbols  have  been  retained  for  future 
reference)  so  that  a  crude  upper  bound  can  be  placed  on  the  phase  bias 
contributed  by  this  term: 


|4>bl  “o|s-s|T  (5-6) 

where  denotes  the  phase  bias  and  T  is  the  signal  duration. 

It  is  desired  to  obtain  an  idea  of  the  typical  size  of  the  estima¬ 
tion  errors  for  the  time-stretch.  The  estimation  errors  for  matched 
filter  receivers  have  been  considered  by  VanTrees^O  and  by  Cook  and 
Bernfeld^^,  who  use  a  Cramer-Rao  lower  bound  approach  to  derive  an 
expression  for  the  minimum  error  variances  for  this  estimate,  and  who 
show  that  the  estimation  errors  approach  the  Cramer-Rao  bound  for 
maximum-likelihood  estimators. 

In  general,  the  estimation  error  for  s  will  be  coupled  with  that  for 
T,  but  a  sufficient  condition  under  which  these  estimates  become 
uncoupled  is  that  the  signal  employed  has  a  real  envelope ^0,  which  has 


already  been  assumed. 
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Adapting  the  discussion  in  Cook  and  Bernfeld,  and  noting  equation 
(2-2),  the  variance  for  the  time-stretch  estimate  is  bounded  as  follows: 


,  (1/m 

E{(s-s)^}  >  - (5-7) 

C  h 


where  h  is  the  input  SNR,  and  C  is  the  RMS  signal  duration,  given  by 


2 


i  f  t2|f(t)|2dt 

^  -T/2 


where  E  is  the  signal  energy,  and  f(t)  is  the  unsampled  signal  envelope. 
Following  the  example  given  in  the  paper  by  Ricker,^  if  the  real 
envelope  has  constant  amplitude,  w^,  and  bandwidth  B>>  y»  the  energy 
becomes 


E  = 


1  2 
2  o 


T 


and 


5 


2 


,  T/2 

-i-  > 

w  T  -T/2 


2 

w 

o 


t^dt 


so  that 


E{(s-s)^}  >  - — -3-- 

(m  T)  h 
o 


(5-8) 


If  the  standard  deviation  corresponding  to  (5-8)  is  considered  a 
useful  estimate  of  the  magnitude  of  the  error  in  the  time-stretch 
estimate,  then  (5-4)  becomes 
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At  -  At(1  ±  — - )  •>  At[1  ±  — - ) 

/h  0)  sT  /h  0)  T 

o  o 

so  that  the  bias  term  is  seen  to  be  negligible  for  most  situations,  where 
/h  u>^T  »  /6.  Using  the  standard  deviation  in  (5-6)  yields 


(5-9) 


This  expression  appears  somewhat  disappointing,  but  in  practice, 

this  large  bias  estimate  is  often  due  to  the  crudeness  of  the  bound, 

rather  than  to  poor  performance  of  the  estimator.  In  many  situations, 

the  time  difference  of  arrival  is  known  to  have  an  upper  bound  (i.e., 

At<At  ),  and  for  a  given  signal,  one  can  use  this  "worst  case"  time 
max 

difference  in  (5-5)  along  with  the  estimate  of  |s-s|  in  (5-8)  to  calcu¬ 
late  an  estimate  of  the  bias,  Indeed,  this  estimator  is  inherently 

limited  to  estimating  delays  that  fall  within  a  range  corresponding  to 
the  time  required  for  the  carrier  to  complete  one  cycle,  i.e. 


or 


min 


+  2ti 

min 


u)  sAt  .  <  0)  sAt  <  0)  sAt  .  +  Zit 

o  min  o  o  min 


where 


At  .  <  At  <  At 

min  max 


6 

.  rain  min 

At  .  =  -  -  - 

min  us  (1) 

o  o 


At  =  At  .  +  — 

max  min  f 

o 
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In  essence,  a  priori  knowledge  of  the  range  of  possible  time  dif¬ 
ferences  is  required  so  that  one  can  assign  the  proper  time-delay  value 


to  an  observed  phase  difference.  If,  for  example.  At  is  known  to  be 

small  and  positive  then  At  .  will  be  selected  as  zero.  At  =  -4— 

min  max  f 

o 

and  the  phase,  i|»,  is  defined  on  (0,2ii). 

If,  for  At  =  At  the  coefficients  f,,  and  f-.  are  nearly  the  same 
max  li  2i 

then  the  imaginary  part  of  the  product  of  the  two  sums  on  the  left-hand 

side  of  (5-5)  will  be  nearly  zero.  For  example,  if  the  carrier  frequency 

is  30kHz,  the  modulation  is  a  850Hz  pure  tone,  the  signal  duration  is 

600  msec,  the  sampling  frequency  is  1700Hz,  and  t  ,  =  0,  At  = 

min  max  i 

o 

then  the  bias  in  phase  is  calculated  numerically  for  real  envelope 
signals  to  yield  the  following  table: 

TABLE  5-1 


Phase  Bias  for  Real  Envelope  Signal 


Input  SNTl 

^fmax  (ysec) 

|(t.bl  (degrees) 

20  dB 

0.0333 

*  0 

0.333 

9.49  x  10"6 

3.33 

1.49  X  10-^ 

33.3  (=  1/fo) 

1.44  X  10”3 

10  dB 

0.0333 

3.43  X 

0.333 

1.55  X  10-^ 

3.33 

1.47  X  10-3 

33.3 

1.47  X  10-2 

0  dB 

0.0333 

1.97  X  10-^ 

0.333 

1.54  X  10-2 

3.33 

1.63  X  10-2 

1 

33.3 

1.65  X  10-1 
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One  would  expect  that  for  more  complicated  real-envelope  waveforms 
with  modulation  frequencies  no  greater  than  850Hz  the  bias  values  would 
be  somewhat  less  than  those  of  Table  5-1 ,  since  for  lower  modulation 
frequencies  the  coefficients  fj^^  and  would  be  more  nearly  the  same. 
Hence,  the  simplified  numerical  computations  yielding  Table  5-1  would 
provide  a  means  of  estimating  the  bias  error  of  receivers  utilizing 
more  complicated  real-envelope  waveforms  of  a  given  bandwidth. 

Because  the  method  used  to  generate  Table  5-1  appears  to  be  of  a 
more  general  usefulness,  the  process  is  parameterized,  and  the  FORTRAN 
source  code  is  included  in  Appendix  C.  The  parameterization  is  useful 
because  with  it  one  can  perform  the  computations  once  for  a  whole  class 
of  signals,  rather  than  having  to  repeat  the  computations  for  each 
specific  signal  used.  This  parameterization  Is  explained  in  the  comments 
of  the  FORTRAN  source  code. 

If  the  signal  used  is  not  a  real-envelope  signal  as  was  assumed 
above  then  the  above  procedure  Is  not  valid,  but  the  bias  can  still  be 
evaluated  numerically  in  the  same  fashion  for  each  specific  signal  using 
the  maximum  time  difference  and  the  standard  deviation  obtained  from  the 
Cramer-Rao  bound.  It  is  mentioned  here  that  If  the  signal  does  not  have 
a  real  envelope,  then  in  general  the  estimates  of  t  and  s  will  in  general 
be  coupled.  In  this  case  the  bound  given  by  (5-7)  will  still  be  valid, 
but  will  not  be  as  tight  as  possible.  A  tighter  bound  that  accounts  for 
T-s  coupling  is  discussed  in  Van  Trees  and  in  Cook  and  Bernfeld. 

It  is  desired  to  combine  the  results  of  this  chapter  with  those  of 
previous  chapters  in  order  to  come  up  with  a  rule  of  thumb  measure  of  the 
total  performance  of  the  estimator.  Here  the  noise  will  be  assumed 
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uncorrelated  (as  in  Chapter  3).  The  measure  of  performance  to  be 

derived  will  be  },  which  represents  the  expectation  of  the 

squared  deviations  of  the  estimator  from  the  true  value,  4)^.  The 

previous  discussion  in  this  chapter  implies  that  the  effect  of  the  t-s 

mismatch  can  be  taken  into  account  by  replacing  with  where 

is  a  zero-mean  random  process  with  variance  determined  from  the 

Cramer-Rao  bound  as  discussed  above.  Without  loss  of  generality, 

will  be  assumed  zero,  so  that  the  variance  of  Figure  3-3  will  now 

2 

represent  E{(<Ji-i()^)  }.  The  densities  of  Figures  3-1  and  3-2  will  here  be 
represented  as  f(4>|4>jj)  In  accordance  with  the  discussion  following 
Equation  3-10.  Denoting  the  marginal  densities  for  the  bias,  i|)^  and  for 
the  total  estimate,  4i ,  as  fC't'^)  and  f(ij>),  respectively,  the  density  f(ij>) 
can  be  expressed  as^ 

Tt 

f(^)  =  /  fC^kb^^^^b^'^^b- 

-n 

Then 

E{((>-<>^)'‘}  =  E{(^^}  =  f  f(<P)d<^ 

-n 

IT  ir 

=  /  {  /  ^  f(^|^j^)d^}f(4.j^) 

-TT  -71 
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where  the  order  of  integration  has  been  interchanged.  Noting  that 
2  2  2 

iji  =  the  inner  integral  becomes 

I  4)  f(<H<t^)d<^  =  /  (4>-4>  )  f(^|4>^)d4. 

-TT  — TT 

2  ^  TT 

-  /  f(^|<|ij^)d<^  +  2^^  /  4>f(4>|<^^^)d^  . 

-1T  -TT 

In  this  equation,  it  is  noted  that  the  first  integral  is  precisely  the 
2 

variance  a  ,  plotted  in  Figure  3~3,  the  second  integral  is  equal  to 

unity,  and  the  third  integral  is  precisely  it>.  ,  so  that 

b 

+  4>^]f(<}>j^)d(t.^^ 

-•TT 


=  .  (5-10) 

Hence,  under  the  assumptions  stated,  one  may  merely  add  the 

2  '  2 
variance  o,  obtained  from  the  Cramer-Rao  bound  to  the  variance  o  from 
D 

Figure  3-3  to  obtain  a  rule-of-thumb  value  for  the  deviations  of  the 
estimator  about  the  true  value.  Note  that  while  this  rule-of-thumb 
was  derived  for  the  uncorrelated  noise  case,  it  is  also  valid  for 
the  correlated  noise  case  if  the  signal-to-noise  ratio  is  high  enough 
that  the  mean  of  the  estimator  is  nearly  equal  to  the  true  value,  <t>  . 


59 


This  concludes  the  discussion  of  the  effects  on  the  estimator  of  a 
T-s  mismatch.  To  summarize,  if  the  time  difference  of  arrival  is  small, 
tiien  the  only  effect  of  a  t-s  mismatch  is  a  slight  bias  as  in  (5-4)  and 
a  decrease  in  the  effective  signal-to-noise  ratio.  If  the  time 
difference  of  arrival  is  not  small,  then  the  phase  becomes  biased.  A 
crude  estimate  of  this  bias  is  given  in  (5-9).  This  estimate  is 
independent  of  the  properties  of  the  signal.  A  better  estimate  can  be 
obtained  by  assuming  a  maximum  time  difference,  and  by  considering  the 
properties  of  the  signal  employed.  For  real  envelope  signals,  this 
process  can  be  parameterized  and  computed  via  the  program  in  Appendix  C. 
For  more  general  signals,  one  must  use  the  specific  signal  to  calculate 
the  phase  of  (5-5).  A  mismatch  in  t  still  affects  only  the  effective 

sigal-to-noise  ratio. 

2 

The  variance,  cTj^,  can  be  added  to  the  variances  given  in  the  graphs 
in  the  previous  chapters  to  provide  a  rule-of-thumb  measure  of  the 
estimator  performance. 
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CHAPTER  6 

SUMMARY  AND  CONCLUSIONS 

It  has  been  shown  in  the  previous  chapters  that  wideband  signals  may 
be  used  for  phase-comparison  time  delay  estimation  provided  the  signals 
employed  have  a  symmetric  power  spectrum.  Chapter  3  analyzed  the 
performance  of  the  estimator  operating  in  uncorrelated  noise,  while 
Chapter  4  discussed  the  performance  in  correlated  noise.  The 
uncorrelated  assumption  is  often  used  in  practice,  where  correlation 
properties  of  the  noise  may  not  be  available. 

It  was  shown  that  in  correlated  noise  the  estimator  becomes  biased, 
but  approaches  an  unbiased  estimator  for  favorable  signal-to-nolse 
ratios.  For  signal-to-noise  ratios  in  which  the  densities  derived  are 
valid,  the  variances  fall  off  as  the  sample  means  and  variances 

approach  tha  derived  theoretical  values.  By  expressing  the  number  of 
input  samples,  N,  as  the  time-bandwidth  product  of  the  signal,  the 
estimator  was  shown  to  use  the  processing  gain  of  large  time- bandwidth 
product  signals  to  reduce  the  estimator  variance. 

In  Chapter  5,  the  effects  of  a  t-s  mismatch  were  examined.  Central 
to  this  is  the  Cramer-Rao  lower  bound  on  the  time-stretch  estimate.  A 
rule-of-thumb  was  given  for  estimating  the  squared  deviations  of  the 
estimator  from  the  true  value. 

This  method  of  time  delay  estimation  has  advantages  over  existing 
methods  in  that  it  can  handle  Doppler  shifted  channels  with  ease,  and 
that  it  can  Identify  and  estimate  the  time  delay  for  resolvable  scat- 


terers  each  w^-h  distinct  t-s  values. 
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further  work  remains  in  evaluating  the  estimator  for  multiple  point 
channels,  and  investigating  possible  interference  phenomena  between  the 
different  point  reflectors.  Also,  work  remains  in  designing  signals  with 
desirable  characteristics  under  the  new  constraint  that  the  signal's 


spectrum  be  symmetric. 
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Appendix  A  -  Correlated  Noise  Model 


This  discussion  of  the  noise  correlation  coefficients  is  taken 
directly  from  Merchant^.  The  noise  model  assumed  here  consists  of  two 
components:  an  incoherent  component  plus  a  coherent  component,  i.e., 


"ki  "ki 


k  =  1,2 


(B-1) 


where  the  superscript  "i"  denotes  the  incoherent  component,  and  the 
superscript  "c"  denotes  the  coherent  component  of  the  noise.  The  vari¬ 
ances  of  the  real  and  imaginary  parts  of  the  incoherent  and  coherent 
2  2 

components  are  and  ,  respectively,  where 

2  2^2 
o  =  o,  +  o  . 
n  i  c 


The  incoherent  noise  component  has  the  following  covariance  properties 


‘'■'ll  ’‘21*  '  ^'>'11  >"21’  ■  °1  °1^  'ij 


E{xj‘  -  Elxj}  yj^)  -  0- 


where  p,  is  the  in-phase  correlation  coefficient  of  the  Incoherent 


noise . 


The  coherent  noise  is  assumed  to  be  the  same  in  each  channel  ex¬ 


c  c 
^21  "li  ® 


cept  for  a  phase  factor: 


If  the  ratio  of  the  coherent  noise  power  to  the  total  noise  power  is 


then  the  covariance  matrix  for  the  Gaussian  density  f(Xj^  ,ypX2,y2) 
becomes 


so  that  the  noise  correlation  coefficients  used  in  Chapter  4  become 

E{x^X2} 

p  - X —  =  r  cos<{>  +  (l-r)p. 

o  2  c  i 

n 


ECx  y  } 

X  =  - X —  =  r  sini^  . 

2  c 

a 

n 


If  the  incoherent  noise  is  uncorrelated  between  channels  (l.e.,  p^  = 
p  and  X  can  be  considered  as  the  real  and  imaginary  components  of  a 

••  ^  r 

"complex  correlation  coefficient,"  y  =  re 
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APPENDIX  B 

Change  in  Correlation  Properties 
Due  to  Complex  Multiplication 


Consider  two  circularly  symmetric  complex  process 


"l 


"2  °  ^2  ^^2 


with  correlations 


2 

E{x^X2}  =  Efy^y^}  =  po 

2 

^^^1^2^  =  -ECx^yj}  =  Xo  . 

If  these  processes  are  each  multiplied  by  complex  numbers  of  unit 
norm,  it  is  desired  to  find  the  new  correlation  properties.  Let 

u^  +  jv^  =  (a  +  jb)(x^  +  jy^)  =  ax^  -  by^  +  jCbx^  +  ay^) 

U2  +  =  (c  +  jd)(x2  +  3X2)  =  CX2  -  dy2  +  j(dX2  +  cy^) 

where  (a  +  jb)  and  (c  +  jd)  represent  arbitrary  complex  numbers  of  unit 

norm.  Then 

E{u^U2}  =  E{(aXj  -  byj)(cx2  -  dy2)} 

=  ECacx^x^  +  bdy^y^  -  adx^y^  -  bcy^x^} 

=  [(ac  +  bd)p  +  (be  -  ad)X]a^  , 
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Similarly , 

E{VjV2}  =  [(ac  +  bd)p  +  (be  -  ad)X]a^ 

ECu^v^}  =  [(ad  -  bc)p  +  (ac  +  bd)X]CT^ 

E{v^U2}  =  ~  ((ad  -  bc)p  +  (ac  +  bd)X]a^ 

E{uiVi}  =  E{u2V2}  =0 

Two  special  cases  are  worthy  of  consideration.  First,  if  both 
processes  n^  and  n2  are  multiplied  by  the  same  number,  then  a  =  c  and 
b  =  d,  so  that 

E{u^U2}  =  E{Vj^V2}  =  (a^  +  b^)po^  =  pa^ 

=  -E{Vj^U2}  =  (a^  +  b^)Xa^  =  Xo^ 

2  2  2  2 

where  a  +b  =c  +d  =1  via  the  assumption  of  unit  norm  multipliers. 
Hence  if  each  process  is  multiplied  by  the  same  complex  number,  the 
correlation  properties  do  not  change. 

Another  special  case  is  that  in  which  only  one  channel  is  multiplied 
by  a  complex  number,  e.g. ,  a+jb=l(a=l,  b=  0).  In  this  case 

E{u^U2}  =  E{v^V2}  =  (cp  -  dX)a^ 

~  E{v^U2}  =  (dp  +  cX)o^  . 

It  is  noted  that  the  general  case  can  always  be  broken  down  into 
consecutive  application  of  the  two  special  cases  (i.e.,  let  c  +  jd  = 

(a  +  jb)(e  +  jf)). 


o  n 
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APPENDIX  C 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


Fortran  Source  Code  for  Phase  Bias  of  Real-Envelope  Signals 

This  program  calculates  the  phase  bias  for  real  envelope 
signals  as  per  pp.  54-55  of  this  thesis  (Hatlestad,  J.D., 
"Phase  Comparison  Time  Delay  Estimation  Using  Wideband 
Signals"  The  Pennsylvania  State  University,  December, 
1985)  . 

The  input  parameters  are  as  follows: 

"time  bandwidth  product"  (FBTEND)  —  This  parameter  is 
the  time  bandwidth  product  of  the  signal.  Since 
the  signal  is  assumed  to  have  a  real  envelope  the 
bandwidth  used  to  calculate  this  must  be  equal  to 
TWICE  the  maximum  modulation  frequency. 

"ratio  of  bandwidth  to  carrier  frequency"  (BWRAT)  —  See 
a>'Ove  note  for  "time  bandwidth  product." 

"ratio  of  time  difference  to  modulation  period" 

(FBDELT)  —  This  parameter  is  the  ratio  of  the 
maximum  time  difference  to  be  estimated  to  the 
modulation  period  of  the  signal  (1/fmax). 

See  page  55  of  thesis. 

"ratio  of  sampling  frequency  to  bandwidth"  (FSRAT)  — 
Must  be  greater  than  or  equal  to  one. 


'input  signal  to  noise  ratio"  (SNR) 
ratio  in  dB. 


Signal  to  noise 


COMPLEX  SUMI , SUMJ,DSUM 
PI  -  4.0*ATAN(1 .0) 

Query  user  for  signal  parameters 
WRITE(6, 10) 

10  F0RMAT($,'  Enter  time-bandwidth  product  of  signal:  ') 
READ(5,15)  FBTEND 

15  F0RMAT(F20. 10) 

FBTEND  =•  FBTEND/2.0 
WRITE(6,20) 

20  F0RMAT($,'  Enter  ratio  of  bandwidth  to  carrier 

Sfrequency:  ') 

READ(5,15)  BWRAT 
BWRAT  -  BWRAT/2.0 
WRITE( 6 ,30) 

30  F0RMAT($,'  Enter  ratio  of  time  difference  to  modulation 
&  period  ') 

READ(5,15)  FBDELT 
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40 


50 

C 

C 

C 

C 


c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c 


c 

c 


c 

c 


WRITE(6 ,40) 

FORMAT($,'  Enter  ratio  of  sampling  frequency  to 
&  bandwidth:  "  ) 

READ(5,15)  FSRAT 
WRITE(6,50) 

F0RMAT($,''  Enter  input  s  i  gnal- 1  o-n  oi  s  e  ratio  (dB):  '  ) 

READ(5,15)  SNR 

Convert  SNR  from  dB  to  decimal 
SNR  =  10.0**(  SNR/10.0  ) 


Compute  auxiliary  variables  (in  radians) 

NP  =•  NINT(  FBTEND*2.0*FSRAT  ) 

FOTEND  =  FBTEND/BWRAT  !  Product  of  carrier  and 

signal  duration 

FBTINCR  =•  FBTEND/FLOAT(NP)  !  Product  of  bandwidth  and 

time  increment 

FOTINCR  =  FOTEND / FLOAT( NP )  !  Product  of  carrier  and 

ti-e  increment 

DELS  =  SQRT(6.0/SNR)/(2.0*PI*F0TEND)  !  Expected  error 

for  s  estimate 

S  =■  1.0  +  DELS 

Convert  arguments  of  trigonometric  functions  to  radians 
FBTINCR  =•  2 .0*PI*FBTINCR 
FOTINCR  =■  2 .0*PI*F0TINCR 
FBDELT  »  2.0*PI*FBDELT 

’ntialize  sums 

SUMI  -  CMPLX(0. 0,0.0) 

SUMJ  =•  CMPLX(0. 0,0.0) 


^ompute  first  sum 
DO  I=‘0,NP 

C0SC0SI=C0S (  I*FBTINCR  )*C0S(  S*I*FBTINCR  ) 

SUMI=»SUMI  +  C0SC0SI*CEXP(  CMPLX (  0 . 0  ,  DELS *I *F0TINCR)  ) 
ENDDO 


Compute  second  sum 
DO  J=*0,NP 

C0SC0SJ=C0S(  J*FBTINCR  )*C0S(  S* ( J*FBTINCR+FBDELT)  ) 
SUMJ=SUMJ+C0SC0SJ*CEXP(  CMPLX ( 0 . 0 , -DELS *J*F0TINCR)  ) 
ENDDO 

Multiply  sums 

DSUM  =•  SUMI*SUMJ 


C 

C  Compute  phase  bias 

PHIB  =  ATAN(  AIMAG(DSUM) /REAL(DSUM)  ) 
PHIS  =  PHI3*130.0/PI 
WRITE(6,60)  PHIB 

60  FORMATC'  PHIB  (DEGREES)  -  ',E10.4) 

C 

STOP 

END 
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